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Abstract 

A toolbox for the development and reduction of the dynamical models of nonequi- 
librium systems is presented. The main components of this toolbox are: Legendre 
integrators, dynamical postprocessing, and thermodynamic projector. 

Thermodynamic projector is the tool to transform almost arbitrary anzatz to a 
thermodynamically consistent model, the postprocessing is the cheapest way to im- 
prove the solution, obtained by the Legendre integrators. Legendre Integrators give 
the opportunity to solve linear equations instead of nonlinear ones for quasiequilib- 
rium (MaxEnt) approximations. 

The essentially new element of this toolbox, the method of thermodynamic pro- 
jector, is demonstrated on application to FENE-P model of polymer kinetic theory. 
The multy-peak model of polymer dynamics is developed. The simplest example, 
discussed in details, is the two peaks model for Gaussian manifold instability in 
polymer dynamics. This type of models opens a way to create the computational 
models for the "molecular individualism" . 
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1 Introduction 



There are many attempts to fill the gap between the microscopic and the macroscopic 
models (the famous micro-macro gap), and to construct closed macroscopic equations. 
Most of the closure assumptions have a relatively narrow domain of applicability, and 
their usage has the following problems: 

1) Violation of the basic physics (thermodynamics) laws; 

2) Absence of the accuracy control procedures; 

3) Absence of the successive step-by-step procedure of the refinement of a model. 
The main object of investigation is the evolution equation 

^ = Jin (1) 

where J is some operator, and \1/ is the distribution function over the phase space. 

The constructed methods are aimed at extracting the dynamics of the macroscopic 
variables from the microscopic equations (^. The prototypes of these methods are the 
quasiequilibrium approximation, dual integrators and the thermodynamic projector. 

The quasiequilibrium closure for the set of the macroscopic variables M (\1/) is built 
with the help of the solution to the variation problem (MaxEnt approximation)^: 

S'(\E') — s> max 

(2) 

M(^) = M, 

where 5'(\E') is the entropy. 

The quasiequilibrium closure is always thermodynamically consistent, but the problem 
2 (the absence of the accuracy control) remains unsolved, and the problem 3 (the absence 
of the refinement procedures) can be solved by adding new macroscopic variables to 
the problem Q. But uncontrolled enlargement of the macroscopic variables set give 
us no guarantee of the accuracy improvement. There exists one more specific problem 

^From time to time it is discussed in the literature, who was the first to introduce the quasiequihbrium 
approximations, and how to interpret them. At least a part of the discussion is due to a different role 
the quasiequilibrium plays in the entropy-conserving and the dissipative dynamics. The very first use of 
the entropy maximization dates back to the classical work of G. W. Gibbs PP, but it was first claimed 
for a principle by E. T. Jaynes j^]- Probably the first explicit and systematic use of quasiequilibria 
to derive dissipation from entropy-conserving systems is due to the works of D. N. Zubarev. Recent 
detailed exposition is given in For dissipative systems, the use of the quasiequilibrium to reduce 
description can be traced to the works of H. Grad on the Boltzmann equation 0). The viewpoint of two 
of the present authors (ANG and IVK) was influenced by the papers by L. I. Rozonoer and co-workers, 
in particular, [Sj • A detailed exposition of the quasiequilibrium approximation for Markov chains 
is given in the book [H] (Chapter 3, Quasiequilibrium and entropy maximum, pp. 92-122), and for the 
BBGKY hierarchy in the paper j^. We have applied maximum entropy principle to the description the 
universal dependence the 3-particle distribution function on the 2-particle distribution function F2 in 
classical systems with pair interactions JHl ■ A very general discussion of the maximum entropy principle 
with applications to dissipative kinetics is given in the review The methods for corrections the 

quasiequilibrium approximations are developed in [171 1381 1121 113| . 
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for the quasiequilibrium approximation Usually while solving the variation problem 
Q we can find explicit dependencies ^'(A) and M(A), where A are the corresponding 
Lagrange multipliers (dual variables), more or less easily. Much more difficult is to find 
the dependencies A(M) and \E'(M) which we need for the closure of the macroscopic 
equations. 

The method of the Legendre integrators consists of building and solving the equations 
of motion for the dual variables. The methods of the first order, based on this idea were 
suggested and tested in the papers O [121 HE] 

The method of the thermodynamic projector let us to represent every ansatz-manifold 
as the solution to the variation problem with the specially chosen constraints. 

The thermodynamic projector is the unique operator which transforms the arbitrary 
vector field equipped with the given Lyapunov function into a vector field with the same 
Lyapunov function (and also this happens on any manifold which is not tangent to the 
level of the Lyapunov function). 

Equations which are derived by the method of the thermodynamic projector are always 
thermodynamically consistent. Although this idea was published in the year 1992 \l7\ . 
the full construction is published only recently in application to the chemical kinetics [TH] . 

One of the problems, discussed in this paper, is to construct the method of the thermo- 
dynamic projector for the derivation of the physically consistent macroscopic equations 
for the polymer dynamics. In the process of building the thermodynamic projector and 
the quasiequilibrium approximation is involved the Lyapunov function for the equations 

which is the entropy S. The equations for the polymer dynamics (Fokker-Planck 
equation) allows us to use the huge amount of different Lyapunov functions and each of 
them can be formally chosen to describe the macroscopic processes. We need to analyze 
the different Lyapunov functions for the Fokker-Planck equation. 

The problem of the accuracy estimation of the resulting approximations and their 
further improvement is suggested to solve with the procedures of the post-processing. 

Suppose that for the dynamical system (0) the approximate invariant manifold has 
been constructed and the approximate slow motion equations \1'm(^) have been derived: 

^ = P^M{^m)), (3) 

where Pi^^,^ is the corresponding projector onto the tangent space T^^^^ of "^m- Suppose 
that we have solved the system Q and have obtained \E'j\/(t). Let's consider the following 
two questions: 

• How well this solution approximates the true solution \l/(t) given the same initial 
conditions? 

• How is it possible to use the solution "^Mit) for it's refinement without solving the 
system Q again? 

These two questions are interconnected. The first question states the problem of the 
accuracy estimation. The second one states the problem of post-processing. 
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The corresponding methods to answer these questions are developed and described in 
this work. 

2 Elimination of fast variables with the help of the 
Lyapunov function 

The most popular way to investigate the dynamics of complicated systems is to split 
the motion into the slow and the fast components, and then to exclude the fast component. 
As a result, one gets a system of equations that describes the evolution of the slow 
variables. The necessary conditions of usefulness of this method arc usually formulated 
as a set of restrictions for the possible dynamics of the "fast subsystem" . Here the "fast 
subsystem" is the subsystem which describes the evolution of the fast variables with an 
assumption that slow variables are constant. 

Unfortunately, often appear situations where wc cannot avoid using this method, and 
there is no proof that it is valid. These situations appear almost everywhere in physical 
kinetics. Here one follows with the same scheme: the relaxation processes are splitted 
into slow and fast. In spite of the fact that in most cases the proofs of vahdity of this 
scheme are absent, the experience helps to avoid fatal errors. 

In this section the method to obtain the equations of the macrokinetics from the 
microdescription is demonstrated. The basis of the analysis is the assumption that if the 
macroscopic variables are chosen in the proper way, then all other variables relax fast: 
the probability distribution of the microscopic variables after a small period of time is 
determined with good accuracy by the macroscopic variables. Let us call this assumption 
the "quasiequilibrium hypothesis". 

The notion "macroscopic variables" is a somewhat relative and is introduced to stress 
the difference of these variables from "everything else". For example, one-particle distri- 
bution function can be "macroscopic" for the full description of the system. 

The goal of this section is to describe the most primitive procedure of derivation of 
the equations for the slow variables and to discuss the form of these equations. 

In this paper the reduction of description goes on with the help of the Lyapunov 
functions. This formalism is the case of the known principle of the conditional maximum 
of entropy with given values of the macroscopic variables. 

Let us review the basic notions of the convex analysis which are used here. 

The subset U of the vector space E is convex, if for every two points Xi,X2 G U it 
contains the segment between xi and X2: for every A e [0, 1] 

Xxi + (1 - X)x2 e U. (4) 

The intersection of any number of the convex sets is convex. 

The convex envelope of the subset M of a vector space E is the smallest convex set 
coM C E, that includes M. It is the intersection of all the convex sets, that include M. 
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If the set U (Z E is convex , xi, Xk G U, Ai, > 0, J2i -^j = then AjXj G t/. 
It leads to another definition of the convex envelope: 

k 

coM = XjXilxi, ...,Xk e M, Ai,...,Afc>0, ^ A^ = 1, A; < 00}. (5) 

1=1 i 

If dimE = n, then in the equation (jS)) it is sufficient to take k < n + 1 (Carthedory 
theorem). 

The function /, defined on the convex set U G E, is convex, if its epigraph, i.e. the 
set of pairs 

Eptf = {ix,g)\xeU,g>fix)}, (6) 

is the convex set in E x R. Sometimes it is convenient to consider functions which can 
reach the value f = 00. If there occurs a necessity to study the functions / which 
are defined on the non-convex set V G E, then it is supposed that / is convex, if the 
restriction of / onto every convex subset of V is convex. If the restriction of / onto every 
line segment from the region of definition is convex, then / is convex. The differentiable 
function / of the class is convex if and only if the matrix of the second derivatives 
d'^f/dxidxj is nonnegative defined (i.e. all its eigenvalues are nonnegative). The smooth 
convex function / on the convex set U C -R" satisfies the inequality 

fix') - fix') > iVfU.,x' - x') = J2idf/dxi),=A^l - x',), ix\ x'eU). (7) 

i 

Geometrically it means, that the graph of / is located above the hyperplane, tangent at 
the point x = x'^. 

The function / is called strictly convex if in the domain of the definition there is no line 
segment, on which it is constant and finite ifix) = const 7^ 00). The sufficient condition 
for the differentiable function / of the class to be strictly convex is that the matrix of 
the second derivatives d'^ f /dxidxj is positive defined. 

In the set of the maximum points of the convex function / on the compact set U iU 
may be not convex) there are some boundary points of U, and if U is convex, then there 
are some extreme points of U. The set of the minimum points of / on the convex set U 
is convex (but may be empty). The strictly convex continuous function has its maximum 
only in the boundary points of U, and if U is convex, then in the extreme points. The 
strictly convex function may have the finite minimum only in one point. 

The function / called concave if the function —f is convex. 

Every bounded convex function on the open subset of R" is continuous. 

Let the C^-smooth function H be defined in the domain U C -R". Let us correspond 
the vector fi = V^i? : /ij = dH/dxt to every point x E U. If the matrix dfii/dxj = 
d'^H / dxidxj is non-degenerated, then for the transform x fi there locally (in the 
neighborhood of every point) exist the differentiable inverse transform. The variables 
fi are often called conjugated variables, and the transform x — >■ /i is called "transition 
to the conjugated coordinates". Let the transform x — > /x be invertible on the open set 
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V (Z U. This means that the function is defined on V. Assuming the smoothness 
of this function, we describe the inverse transform /x — > a; in the same way as the direct. 
For this purpose we introduce a function 



The function G called the Legendre transform of H. 

With the help of the conjugated coordinates it is possible to write down the necessary 
conditions of the extremum for the problems with the linear constraints on the open set 
in a very simple way: 



With the method of Lagrange multipliers we get the system of the equations which is 
giving us the necessary conditions for the solution to the problem (0): 



where the A, are the Lagrange multipliers. The necessary conditions for the extremum 
are given by the system of the equations (fTUI) . One part of the system is linear in the x 
coordinates, and the other part is linear in the conjugated coordinates /i. 

Let us have the Legendre transform for the function H{x), let the transform 

X ^ fi have the smooth reverse, and let the solution to the problem Q be unique for 
some open set of the values of the vector {Mi, Mk) G -R". Also let the point of the 
minimum Xmin, and, consequently, the minimal value of H be smooth dependent on M, 
Hrain = H{M). Let us denote /im, = dH{M)/dMi, hm = (/iMi, /^mJ- Let us get 
some information about the function H{M) from the functions H{x) and G{n) without 
solving any equations. With the known value of the vector fiM we can immediately find 
the vector /i in the corresponding point of the conditional minimum, fij = J2i f^Mi^ij. 
From this equality we get 



G{fi) = {fi,x{fi)) - H{x{fi)) = '^fiiXi{fi) - H{x{fi)), 
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H{x) min, 

X],- rriijXj = Mi, {i = 1, ...,k), X e U. 



(9) 



fij = Y.. XiUiij, j = 1, ...,n, 
Y,. rUijXj = Mi, (i = 1, k). 



(10) 



xifiAi) = (V^G(/i))|;,^=j:,^^,.m,,. 
From Xf^j^j we obtain M{^m) and H{M{jj,M))'- 

Mi{nM) = ^mijXj{^iM), H{M{^m)) = H{x{fiM))- 



(11) 



(12) 



j 



Finally, the Legendre transform G{^m) for the function H{M) is: 



(13) 
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So, we can find dependencies fi{^M), x{^m), M{^m), H{^m) and G^hm) from tlie func- 
tions H{x) and G{x) without solving any equations. We hope, that the similar notations 
for H{x) and corresponding conditional minimum function H{M), and for their Legen- 
dre transforms G{fi) and G{^m) will not cause a confusion. Let us note, that with our 
assumptions the reversibility of the transform M fiM follows from the reversibility of 
the transform x ^ fi, and moreover, the function M{^m) can be found explicitly. 

The convexity of the function H{x) usually makes our assumptions (existence and 
uniqueness of the conditional minimum, global reversibility of the transform a; — > /i, 
smoothness of the function H{M)) easier to check. Note, that the convexity of the 
function H[M) is neither necessary nor sufficient condition for our assumptions. If H{x) 
is convex, then the function of the conditional minimum H{M) is convex too. 

Now we proceed to the problem of elimination of the fast variables. Let us have the 
system of differential equations with smooth right-hand sides 



in the convex domain U C -R", and moreover let the linear transform x — > M, Mj = 
rriijXj from the phase space to the space of the slow variables M be defined. We can 
assume, that we have no linearly dependent rows in the matrix rriij, because it is always 
possible to eliminate the linear dependent functions Mj(x), if they are present. 

Let us assume that in the interesting for us domain of the initial conditions xq the 
solutions x{t) of the equations p4|) are developing in the following way: the vector x{t) 
is going fast to the value which is defined by the slow variables M; after that x can be 
represented as the function of M with a good accuracy, and this function is unique for 
every initial conditions. So, 

A) For each value of the slow variables M G M{U) there exist such x = x*{M), that 
if M{x^) = M°, then x{t) is going very fast to some small neighborhood of the x*(M°), 
and during that M(x(t)) is almost constant. 

B) In the process of the further evolution, x{t) stays in the small neighborhood of the 
value of X which corresponds to M{x(t)), so x is close to x*{M{x(t))). 

It is usually impossible to give a strong proof for A and B for the situations of real 
complexity in the nonequilibrium thermodynamics, so this assumptions are, probably, the 
weakest point of all the construction. We are accepting them because we are sure that the 
evolution of the macroscopic variables is possible to describe by the autonomous system 
of differential equations of the first order (if it is impossible, then, probably, one should 
extend the list of the macroscopic variables with respect to the physical properties of 
the investigated process). There is another way to deal with this problem: to equip the 
obtained approximations by the postprocessing. The postprocessing helps us to correct 
the errors, if they are not too big, and gives us a signal if they are too big. 

If we know the function x*{M), then we can write 



X = F{x) 



(14) 



M = mF{x*{M)), Mi = "^nHjFjix* {M)). 



(15) 



j 
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In general, this equation can be used only for short periods of time which do not exceed 
some limit. The right-hand side mF{x*{M)) of the equations ()15|) is not exactly the 
mF{x{t)), and it may cause the error increment, and as a result the solution of the 
equations (jl5|) will divert from the true solution strongly. The exclusion is the case when 
in accordance to the equations (|15|) M{t) tends to the only stable fixed point when t — > oo. 
If the solution of the equations (fT^ and the real values of M{x(t)) are not succeed in 
getting far one from another during the time in which the solution of the equations (jlSp 
is coming in the small neighborhood of the fixed point, then the equations ()15|) can be 
used also for t — > oo. 

The function x*{M) for the particular system is not unique, but the range of choice 
is small in that sense, in which the neighborhood of x*{M{x{t))) (in which the evolution 
goes after the short period of time) is small. 

Let us have the Lyapunov function H{x) for the system ()14j) which is decreasing along 
the trajectories. We can try to find the dependence x*{M) as the solution to the problem 
H{x) — >■ min, mx = M. This way seems to be natural, but it does not follow directly 
from the assumptions A and B. For example, there could be a situation in which H is 
very sensitive to small changes of the slow variables, and not sensitive to the changes of 
the fast variables. In this situation the assumption, that x* (M) is the point of conditional 
minimum of the function H, may not give the desired result. The following idea does 
not solve the problem, but it can be useful: In the applications, the system (fT^ usually 
dependens on some parameters. It seems to be more reasonable to use the Lyapunov 
function which does not depend on these parameters, if there exists such a function. It is 
most important in the case when among the parameters we have such, that their values 
are determining, whether is it possible to split the variables to fast and slow, or not. 

So, the fast variables will be eliminated with the help of the Lyapunov function. Let 
us have the Lyapunov function H for the initial system, let the transform x ^ fi = VxH 
have the smooth inverse, and let us know the Legendre transform G{fi) for the function 
H{x). Here it is also assumed that for every M G M{U) the problem (jU)) has the unique 
solution, and the minimum point x*{M), and the function of the conditional minimum 
H{M) smoothly depend on M. With the value fiM = ^ mH{M) it is possible to find 
fi{fiM), x{fi{fiM)) (look at the (jlllll3j) ). The result is 

M = mF(V^G(/i))| (16) 

where is the product of the row vector and the matrix m: 

i 

Vf^G is the vector with the components dG/dfii, and all derivatives are taken in the point 
fJ' = fJ'M^- The right-hand sides of (fTB|) are defined as the functions of fiM- In order 
to define them as functions of M, one needs to make the Legendre transform, find the 
function H{M) and, respectively, fiM = ^ mH{M) from the function G{hm) (US))- It is 
impossible to make these calculations explicitly in such a general case. It seems to be 
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a very natural and convenient to define the right-hand sides of the kinetic equation as 
the functions of the conjugated variables. If in the beginning the right-hand sides of the 
equation (fT^ are defined as the functions of /i (i.e. x = J{fJ.)), then the equations (fTH|l 
have a very simple form: 

M = mJ{nMm). (17) 

H[M) is the Lyapunov function for (jl6|) . its time derivative due to the system (jl6|) is not 
positive: 

H{M) = {fiM,mJ{fiMm)) = {fiMm, J{fiMm)) < 0, (18) 

because (/i, J(/i)) = H{x) < 0. 

Let us call the systems dissipative, if if < and conservative, ii H = 0. For the 
dissipative system we have H{M) < (jl8|) . and if the system is conservative, then for 
all values of fi we have (/i, J(/u)) = H{x) = 0, and then from the equation ()18|) we get 
H{M) = {fiM,^J{fJ'M^)) = ifJ'M^, J{i^mi^)) = 0. So, we proved the following 

Theorem^. The Lyapunov function for the microscopic system remains the Lyapunov 
function for the macroscopic system, and if the microscopic system is conservative, then 
its projection to the space of the macroscopic variables remains conservative. 

If necessary, it is easy to perform further exclusion of the variables in the equations 
()16|) with the help of the function H[M). The right-hand sides of the resulting equations 
will be defined again as the functions on the conjugated variables, and the function of 
the conditional minimum will be the Lyapunov function again. Let us note that in (jl7|) 
we have neither H nor G in the explicit form (they occur only when we need to find the 
connections between M and or x and /i). 

Convexity of H was never used above, but the natural domain of applicability of 
the described formalism are systems with convex Lyapunov functions H{x), or at least 
with such H, that the sets {x\H{x) < h} are convex. Otherwise there exist such linear 
manifolds, that the local minimum of H is not unique on them, and further considerations 
are required to select the relevant minima. The finite dimensionality of the phase space is 
not so important, because everything said above can be applied to the infinite-dimension 
case with proper restrictions. Let E be the Banach space, U G E he the convex open 
set, H : U Rhe C^-smooth function. With every point x E U we associate the linear 
functional fi^ G -E*: Hx = ^ xH which is the differential of H in the point x. Let V be 
the set of the values of for x E U and let us have the smooth mapping J from E* to 
E in the neighborhood of V. The system {U, H, J) determines the system of equations: 

X = J(/ix). (19) 

Let L he the closed subset of E and for every M E U/L let the problem H{x) 
min, x/L = M, x E U have the unique solution Xmin which is C^-smooth dependent on 
M, H{M) = H{xmin)- Denoting /im = VmH{M) G {E/L)* G E* we can define the 

^This is a rather old theorem, one of us had published this theorem in 1984 already as textbook 
material ([S], chapter 3 "Quasiequilibrium and entropy maximum", p. 37, see also the paper ^3), but 
from time to time different particular cases of this theorem are continued to be published as new results. 
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factor-system which is the exact analogue of (fT^ : 



M = J{^m)/L. (20) 

Here the argument J is the hnear functional on hm'- I^'MX = fiuix/L) 

The described procedure of the elimination of variables has one very important com- 
mutativity property: If one makes a further simplification and transact to the variables 
= N{M), then after the application of the described formalism to the system (|2U|1 
with the function H{M), one get the same result as after the application of this formal- 
ism directly to the reduction from x to N{x) = N{M{x)). So, the chain of exclusions 
X ^ M ^ N gives us the same result as the direct exclusion x N. 

3 The main problems in usage of the quasiequilib- 
rium approximations 

Our problem is to build the closed system 

M = J{M), 

from the initial system (P) and its Lyapunov function. 

If we know the function x*{M) then it is sufficient to calculate m{F{x*{M))). This 
problem is the problem of the calculation of the projection of the microscopic vector field 
F on the macroscopic variables M in known point x*{M). Let us call this problem the 
problem about the macroscopic projection. If the right-hand parts are expressed through 
fi then we have the problem about the macroscopic projection too. 

Another problem is to find fiM- Usually it is necessary to solve the system of non- 
linear equations (if the function H is not quadratic) to solve this problem. Indeed, let 
us consider the conditions for the conditional extremum of H with given values of the 
moments M. From the functions H{x), G{fi) we get fj,{fiM), x{fiM), M{fiM), H{M{iim))- 
But in this list we have no function ^m{M). We can find this function as the solution of 
the equation 

M{fXMm) = M. (21) 

Let us give a few examples. 

One-particle approximation. Let x be the A^-particle distribution function, /7v(^i, C,n), 
where is vector of coordinates and momenta of the i-th particle, and let the evolution 
of this function be described by the linear equation 

Furthermore, let M be one-particle distribution function 



= Nj fNi^,^2,...,^N)d^2:.d^N, (23) 
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and H be the entropy (we use the if-function which is equal to the minus entropy) 

HUn)= J /iv(ln/iv-l)ci^e, (24) 
For given f^, H, f we get /i = In f^, Jn = exp /i, 

G(/i)= / exp/i(ei,...,e^)ci^e, (25) 



^{Jn) = J J2iLi ^(^ ~ CdfNi^i, ^N)d'^^'-i the extremum conditions (|Tn|l are of the form 

TV 

Mei,-,eTv) = / ^^e/ii(oE'^(^-^^) = E/^^^^^)' 

i=l i 

/jv = exp^/ii(^i). (26) 

i 

The normahzation condition here is J fNd^^ = 1, that is 

j exp/ii(Orfe = l- (27) 

Connection between the macroscopic variables /i (that is M) and the quasiequilibrium 
values of the microscopic variables (that is, x^^) is given by well known formula: 

/Tv(6,-,eTv) = ^/i(6)-/i(eTv). (28) 



Projection of the microscopic vector field ()22j) can be found by direct integration. 

Two-particle distribution function as the macroscopic variable. One-particle distribu- 
tion function fi{C,) is often not sufficient because, for example, the energy of the interaction 
of pairs of particles cannot be found from this function. Much more detailed description 
is given by the two-particle distribution function 



/2(ei,6) = iV(iV- 1) J f^{^,,...,^N)d^,...d^N. (29) 
We can easily find the expression 

Ai(6, ■■■,^n) = ^ /^2(^i,0); 
/7v(6,--,^jv) = exp/i = exp ^ /i2(6,'^j)- 

But it is difficult to find the connection between /i2 and /2 explicitly. Only a series 
expansion for it in the neighborhood of the uncorrelated state is known jTUI. The problem 
about the macroscopic projection becomes hard too: the necessary integrals in general 
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case are impossible to find analytically. For two-particle distribution functions as well as 
for majority of the most interesting variables the transform M ^ /zm is very complicated 
in the direct direction and not very simple (as simple as the derivation of /2 from f^) in 
the opposite direction. 

So, we need to avoid the necessity to calculate hm^M) (and, if possible, to make less 
calculations to find M{ixm))- 

The first of these two problems (avoiding calculation of fiM^M)) is solved by the 
method of the Legendre integrators which is developed by us [H] . 



4 Legendre integrators 

The main idea of the Legendre integrators is to find some alternate way to solve the 
macroscopic equations M = J{x): a way to find their solution in the absence of the 
explicit form of these equations. 

First of all, note, that we have a linear connection between M and fiM'- 

^ = imiDlSix)rV)'j^; (30) 



M = —{■mx{^M'm)) = m{D^x)rnF 11] 

(31) 

D,x = (Z},.^)-i = iDlSix))-\ 

Calculation of the functions m{F{x)) is the standard problem of the macroscopic 
projection. Dependencies x{^m) are usually quite simple. We suggest the following 
advancing in time to solve (unknown) equations M = <I>(M): 

IJM{t) ^x = xifiM) ^M^fiM^ /iAf(t + At) ^ M{t + At). (32) 

In the sequence (|32p there is one operation of macroscopic projection and one operation 
of solving the system of hnear equations (jHU)) . 

Formally, it is possible to write down the equations for fiM- 



dfj, 



M 



dt 



{m{DlS{x))-V)-^mF{x), (33) 



where x = x^^. 

Nevertheless, explicit inversion of the operator in the right-hand part of the equation 
(jHHj) is usually difficult and one should use the chain of computations ((221) • c>^^ ^^st 
calculations using of the Legendre integrators jT^ IT^ the methods of the first order of 
accuracy were used. This is not the principal restriction: the scheme fj32|) gives us a 
possibility to calculate fiM for any given /xm, so all known methods of the higher order 
can be used (for example, the Runge-Kutta method with the different procedures of the 
automatic step selection [221 12E1 12Z|)- 
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5 Lyapunov functions for the Fokker-Planck equa- 
tion 

The Fokker-Planck equation (FPE) in the absence of the drive forces has the form 

^^^^ = V,{D{^{q, t)V,U{q) + V,v[/(g, t))}, (34) 

where \1/ is the probabihty density over the configuration space, q is the point of this 
space, '^{q) is the function of the time t, U{q) is the normahzed potential energy {U = 
U potential I kT), D{q) is positively semidefinite diffusion operator {{yi,Dy) > 0). 
The FPE has two important properties: 

1) Conservation of the total probability: 

^ j v[/(g,t)rfg = 0. (35) 

2) Dissipation: for every convex function of one variable h{a) {h"{a) > 0, a > 0) the 
following functional S[$] is monotonically non- increasing in time: 



dq, (36) 



'(q). 

where 

^*(g) = const ■ exp(-f/(g)), (37) 

is the Boltzmann-Gibbs distribution. 

For h{a) = alna, the functional S[$] is the usual Boltzmann-Gibbs-Shannon entropy: 

sm = -l (38) 

Let us calculate the time derivative of S[^] due to FPE (p^ . Note, that 

^(g) \ _ V,^(g) + ^{q)VgU 



V^*(g) J ^*(g) 
so we can rewrite FPE as follows: 



dt ' \ \^*iq). 

Let us consider FPE in the domain Q. Function dS/dt consists of two summands: The 
first is the integral of the local "production of S" , J cr{q)dq, and the second is the flow 
through the boundary of the domain Q: 



dt Jn V ^* / n V n ^* 

r / * \ / vi/ M 

div 
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where dw is the differential of the area, Uq is a vector of the unitary normal to d^l in the 
point q, a{q) is the entropy S production: 

<^fe) = *-ft"(^)(^.i'V,(^))>0, (39) 

Let the flow of through the boundary dQ be equal to zero: 

Uq, DVq (^11=0 



at all points of dVt. Then 

dS_ 



a{q)dq > 0. 



The most important cases of S selection are: 

h{q) = alna, S" is the Boltzmann-Shannon-Gibbs entropy; 

h{a) = a In ax — a In ax is the maximal family of additive trace-form entropies 1201 
E] (these entropies are additive for composition of independent subsystems); 

= TZ^jCK 7^ 1 is the Tsallis entropy |22]- These entropies are not additive, but 
become additive after nonlinear monotonous transformation. This property can serve as 
definition of the Tsallis entropies in the class of generalized entropies ^T] . 



6 Macroscopic variables and quasiequilibrium distri- 
bution functions for FPE 

The set of the macroscopic variables can be continuous or discrete. Let a be the discrete 
or continuous parameter, that enumerates the macroscopic variables, and Ma be the 
corresponding variables. Every macroscopic value is defined by its microscopic density 
rriaiq) ■ 

Ma= [ ma{q)^{q)dq (40) 

The choice of the domain fi, in which we are solving the FPE, needs to be discussed 
separately. We can suppose formally, that f2 = i?", but for the calculations it is better 
to make it as small as possible with the preservation of the accuracy. Usually, when 
||g|| — >■ oo the function \E'(g) tends to zero faster, then exponential, and we can a priori 
select the bounded domain f2, out of which ^ is negligibly small. 

We shall do the calculations for the general form of 5* (see equation (j36|) ) and give the 
examples for the most popular choice (jHHj) of S. 

Quasiequilibrium function Ma for the given Lyapunov function S (IHUj) is defined as 
the solution to the problem 

/ S{^!)^ max 

\ jma{q)'^{q)dq = Ma ' ^ ' 
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Due to the convexity of h (and, consequently, concavity of S), it is sufficient to inves- 
tigate the conditions of the local extremum: 



^vj-S" = ^m„(g)/i„, (42) 

a 

where /i^ are variables, dual to Ma (fJ'M)- For continuous parameter the sum in the 
equation (jl^ is replaced by integration on a. 

Next, we use the standard Riesz representation of functionals (through the scalar 
product). Let us write 



For the quasiequilibrium distribution we have 

^ = ^*g(^-Y,^o.{q)fi}j, (43) 

where g{a) is a function of one variable, inverse to h'{b). Note, that h'{b) is a monotonous 
increasing function (because h is convex), so g{a) is a monotonous increasing function 
too, and g'{a) = {h"lg{a)))-\ 

Let us denote the quasiequilibrium distribution function P^j) as \E''^'^({/iQ,}, g). 
For the BGS entropy h{b) = b{\nb — 1), h'{b) = In 6, g{a) = expa, and the equations 
transforms into the following equation: 



^''({/Wa},g) = ^*exp ^-5^m„(g)/i„j 



(44) 



For the next steps it is convenient to consider the temperature dependence explicitly 
(i.e. write f3U instead of U in FPE, (3 = 1/kT), then we have \E'* = const ■ exp(— /5[/). 

For the classical BGS entropy the quasiequilibrium distribution will take the 
simplest form: 

^^^({/"a}, q) = exp ^-/io - fiuU - rnaiq)fi^ , (45) 

where = (3 = 1/kT, fiQ is a variable, conjugated to Mq = J^'^dq = 1. The function 
()45|) is a solution to the problem: 

— \I' In "^dq max 
Mo{^) = J^mdq = 1 

Mu{^) = U{q)^{q)dq = Mu ^ ' 

Ma{^) = /^m,(g)^(g)rfg = M„ 



16 



In the problem ()46p we move from the relative (so-called Kullback) entropy to the absolute 
entropy. 

Selection of the macroscopic variables is the most critical point in construction of the 
quasiequilibrium approximations. It is always necessary to select them, basing on the 
specific of the problem. Nevertheless, there are some simple general recommendations 
about construction of the set of variables for the Legendre integrators. 

1) It is necessary to include Mq in the list of variables, because /iq is not constant in 
time; 

2) It is useful to include Mu in the list of variables. With this variable in the process 
of the relaxation all other /Iq, — 0, and iiu — 1/kT. 

3) It is better for the set of the functions ma{q) to be hnearly independent. 
For the classical entropy we have 



^'''^({/i},g,t) =exp -fioit) - fiu{t)Uiq) -J2^c.iq)fia{t) 



Due to the equation (jlTj) we have 



9$ 

= -vp 

dt 



The FPE gives us 
'dt 



(/ic7-/?)(V,DV)f/(g)+ 



^/i„(V,Z^V)m„(g) - ^(2/ic;/i, -/?/!„)( Vf/(g),DVm,(g) ) 

a a 

^ /ia/Xa' ( Vm„ (g) , DVm^, (g) ) 



(47) 



(4J 



(49) 



To calculate ^({/i}) means to calculate the following integrals: 



dMu _ 
dt ~ 
dMa 
dt 



ma{q) dq, 



where ^ is calculated due to equation dM^/dt = 0. 

From the equation (j48|) we get the conditions for derivation of [i 



Mr 



dfi 



u 



-M, 



-Mu 
djj,o 



dt 

dt dt 



dt 

2\ dfJ^u 



dt 



Mn = 0; 



/TT \ dfJ'U 



a 

djjL^ 



dt 



(50) 
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where by {f{q)g{q))^ we denote the averaging {fg)^ = f{q)g{q)-^{q)dq. 

We get the closed system for derivation of the dynamics of /x. But the question about 
the choice of the macroscopic variables still remains open. 

In the problem of the quasiequilibrium we find the projections of \& to the given set 
of the functions (linear space), afterwards we calculate \l/ due to the maximum entropy 
condition. 

It seems to be physically sensible to choose the additional variables to Mq, Mu as the 
projections of \I' onto some equilibrium states: 



There are two classical choices of macroscopic variables: 

1) a = R+ (Laplace transform of the energy distribution density) 

2) a = ik, k E R (Fourier transform of the energy distribution density). 

The variable Mu is the average energy in the potential well U{q). In analogue to 
this, the variable Mq,(\1/) ()51|) for the real a > can be considered as the energy in the 
potential well e~°'^^'^\ This potential is gained by the monotonous nonlinear deformation 
of the energy scale U e-"^^'?). For imaginary a this nonlinear deformation is given by 
the periodical functions U cos{kU) + ism{kU) 

A benefit of usage of ()5H) is also in that {menial) = Mq,+q,', and we have to perform 
less calculations in ()50|). This set of the deformed energies can be used for both the initial 
potential U and the set of additional potentials. 

Is this set of macroscopic variables sufficient for description of nonequilibrium kinet- 
ics of polymers in presence of flow? Probability densities for all the quasiequilibrium 
distributions which can be constructed with this macroscopic variables have the form 
\E'(g) = (p{U{q)), where ^p{U) is a function of one variable. Is this class of distributions 
sufficient for the specific problem? This question can be answered only after specification 
the problem. But what is possible to do, if the closure with these variables gives too big 
error (the estimation of accuracy is discussed below)? There are at least two ways: to 
extend the list of variables or to improve the quasiequilibrium manifold EH] (appli- 
cation of the methods of invariant manifolds to improving the quasiequlibrium closure 
for dynamics of dilute polymeric solution is presented in [H^). The extension of the list 
of variables is the central method of the extended irreversible thermodynamics [lU]- It 
is possible to combine the potential energy U{q), the vector of the configuration space 
q, and the gradient of U{q), VU{q) = —F{q), (-F(g) is the force) and to obtain a huge 
amount of densities m{q) which can be scalars, vectors, or tensors. The corresponding 
"macroscopic variables" are j^m{q)'^{q)dq. 

The best hint for a choice of new macroscopic variables is the analysis of the right 
hand side of dynamic equations [31] . The well known distinguished macroscopic variable 
associated with the polymeric kinetic equations is the polymeric stress tensor |521 
This variable is not the conserved quantity but nevertheless it should be treated as a 
relevant slow variable because it actually contributes to the macroscopic (hydro dynamic) 
equations. Equations for the stress tensor are known as the constitutive equations, and 




(51) 
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the problem of reduced description for the polymeric models consists in deriving such 
equations from the kinetic equation. 
The tensor 

Tp,, = k^T (^6ij - j Fiqj^{q)d(^ (52) 

gives a contribution to stresses caused by the presence of polymer molecules for unit 
density. Here F{q) = —VU{q) is the force vector, Sij is the Kronecker symbol. For 
spherically symmetric potentials {U{q) = u{q^)) this tensor is symmetric. The tensor of 
dencities mij{q) = Fi{q)qj is the first addition to the dencities which depend only of U{q). 



7 Macroscopic variables and boundary conditions 

There is a standard technique to solve the boundary value and initial-boundary value 
problems of mathematical physics: first to build the space of the functions which satisfy 
the boundary conditions, and then to find the solution in this space. 

When one uses the Legendre integrators, a special technique is needed to satisfy the 
boundary conditions. 

FPE describes the evolution of the probability distribution. It conserves the total 
probability. The natural boundary conditions for the FPE is the absence of the flow 
through the boundary of Q: 

vl/^^(z.„DV,(|^)) =0, (53) 

on dQ, where Ug is a vector of outlet normal to dQ in the point q. 

Quasiequilibrium distribution functions (jTTj) satisfy the condition if 

r (z/„DV,f/(g))=0 

\ (z/„DV,m„(g)) = ' ^""^^ 

for all a. 

There is also a different way to satisfy conditions (j53|) : to make '^*\qq = 0. It is 
possible to do by making U{q) — oo while q qo E dQ. But this choice leads to the 
singularities and is very inconvenient from the numerical point of view. 

Conditions (jHl]) look somewhat surprisingly, if considered without the context of the 
quasiequilibrium approximations: for the quasiequilibrium solutions the absence of the 
flow through the barrier follows not from the infinite heights of the barrier, but from the 
fact, that the normal derivatives of U and are zeros. 

To satisfy the condition (jMjl it may be necessary to deform the initial potential U 
and densities m{q). This deformation will be the smoothing of U near dfl. The error, 
introduced by this deformation is usually not very big (because of the smallness of 
near dfl) and can be estimated easily. 
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So, the quasiequilibrium approximation and the Legendre Integrators of any order of 
accuracy are built, and the way to satisfy the boundary conditions is suggested. First 
numerical experiments ^1 1^ proved the effectiveness of this idea. 

The main computational challenge in this method is to calculate the integrals of the 
form 

(Y1 ^kVkiq)^ exp 7i^i(g)) dq (55) 

where (pk{(l),i'i{q) are known functions (usually they are given analytically). For the 
problems of the polymer physics the complexity of the problem fl55|) is dependent on few 
characteristics: 

1) The quantity of the different functions ^Pk{q),4'i{q) is usually 5-10; 

2) The dimension of the space in which the integration goes is usually 10-100. 



/ 

Jn 



8 Thermodynamic projector and Galerkin approxi- 
mations 

Almost every manifold of the functions can be represented as the solution to the quasiequi- 
librium problem (PT|) . if this manifold is not tangent to the level surface of the entropy 
S = const ^\ . For this representation only the right system of restrictions is needed. By 
the simple parameterization with the moments M(\l') it is possible to get only the clas- 
sical quasiequilibrium manifolds The restrictions which are necessary to represent 
manifold Q as the quasiequilibrium manifold are built as follows. Let / e f2, and Tf be 
the tangent space to Q in the point /. On the space of the distribution functions E we 
define the projector Pf : E Tf. Operator Pf depends smoothly on the point / and on 
Tf. The problem of the quasiequilibrium is posed as follows: 

S{-^) max 

Pf(^ - /) = ^^^^ 

The necessary and sufficient condition for / to be the unique solution to the problem (j^Ujl 
is [T?!: 

kei Pf CkeiDSlf, (57) 

that is, if P/(v9) = 0, then DS\f{ip) =0. For the classical entropy /^^[/(v?) = — J ip {q) In f{q)dq 
and the condition (j37|) takes the form: 



If Pf{^) = 0, then j ^In fdq = 0. (5^ 



Among all projectors which satisfy the condition ()57|1 there is unique projector which 
has the following property: let us have the appropriate equation 

^ = Jin 
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for which dS\^]/dt > 0. Then for the projected equation on Q 



f = Pf{Af)), 



(59) 



we also have dS[f]/dt > 0. 

This projector was introduced in the paper [TH], and there it is also proved its unique- 
ness. It is built as follows. 

Let us require that the field of projectors, P(\1/,T), is defined for any \Ef and T, if 



From these conditions it follows immediately that in the equilibrium, P(\I'*,T) is the 
orthogonal projector onto T (orthogonality with respect to entropic scalar product (I)**)- 

The field of projectors is constructed in the neighborhood of the equilibrium based 
on the requirement of maximal smoothness of P as a function of g^jj = D^S and It 
turns out that to the first order in the deviations — vp* and — g^;,*, the projector is 
defined uniquely. Let us first describe the construction of the projector, and next discuss 
its uniqueness. 

Let the subspace T C E, the point and the differential of the entropy in this point, 
= DifiS, be defined such that the transversality condition (pTUjl is satisfied. Let us 

define Tq = Tf^ker^f,^. By the condition (IHUjl . Tq ^ T. Let us denote, eg = eg(T) G T the 

vector in T, such that Cg is orthogonal to Tq, and is normalized by the condition g{eg) = 1. 

Vector eg is defined unambiguously. Projector Ps,^ = P(\1/,T) is defined as follows: For 

any z E E, 



where Pq is the orthogonal projector on Tq (orthogonality with respect to the entropic 
scalar product (|)«(). Entropic projector depends on the point q through the 
dependence of the scalar product (|)^., and also through the differential of S" in the 
functional g^. 

Obviously, P(/) = implies g{f) = 0, that is, the thermodynamicity requirement is 
satisfied. Uniqueness of the thermodynamic projector (|HT|) is supported by the require- 
ment of the maximal smoothness (analyticity) (THj of the projector as a function of gq, 
and (1)^, and is done in two steps which we sketch here: 

1. Considering the expansion of the entropy in the equilibrium up to the quadratic 
terms, one shows that in the equilibrium the thermodynamic projector is the or- 
thogonal projector with respect to the scalar product (|)**. 

2. For a given g, one considers auxiliary dissipative dynamic systems which satisfy 
the condition: For every \E'' G U, it holds, 5',p(J(^')) = 0, that is, g^^, defines 
an additional linear conservation law for the auxiliary systems. For the auxiliary 
systems, the point \1/ is the equilibrium. Eliminating the linear conservation law g^j,, 
and using the result of the previous point, we end up with the formula (|HT|). 



(60) 



PsAf)=Po{z)+eggM 



(61) 
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Thermodynamic projector allows to use almost arbitrary manifolds as quasiequilib- 
rium closure assumption. If the projection of FPE ()59p is built with the thermodynamic 
projector, then dS/dt conserves (not only the sign, but also the value). The only restric- 
tion is that the manifold must not be tangent to the level surfaces of S (and must contain 
the equilibrium point). 

Let us write down the explicit formulas for the closure assumption of the form 

f{q) = ^\q) + Y.Uq)^^c.. (62) 

a 

Due to probability conservation for all a we have / fa{q)dq = 0. 

Tangent spaces to the manifold (jU^ in all points coincide and have the form T = 
{J2a f^afoiiq)}- The natural coordinates in T For every f{q) of the form (j62|) 

there is the entropic scalar product, defined in T: 



/(?) 

In the coordinates /!„ this scalar product has the form 

a p a,f3 

where 

flMIM,, 
= J ^teT 

We will need the orthonormalized basis of the subspace T f]keT{DS\f). This subspace 
is defined by the equation 



J J]/a(g)/ialn^i^f/g = 0. 



Let be J fi{q) In -^r^dq ^ bor the definiteness. Suppose for a > 1 

qa = fa- I^a/l, (63) 



where Ua 



!Uq)^^Mr/Q 

IfM'^^M)dq 



Let us orthogonalize the family of the vectors (a > 1) with respect to the scalar product 
(■|-)/. We will get the orthogonal basis in T f]keT{DS\f): {ea}(a > 1). 

Let ei E The the vector, orthogonal to all Cq, (for example, ei = a,(/i— X]a>i ^aifil^a) /)) 
and let ei be normalized in the following way: J ei{q) In^^dq = 1. The projection of 
the vector J on T is defined in this way: 

PP ^ / Ji,) .„ A|,, , g I M^,,. (64) 
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Projector ()64|) allows to consider every manifold of the form ()62|) which is not tangent 
to the level surface of the entropy S, as the quasiequilibrium manifold. If the vector field 
is projected with the operator then the dissipation is conserved. 

As we can see, there is a "law of the difficulty conservation": for the quasiequilibrium 
with the moment parameterization the manifold is not explicit, and it can be difficult 
to calculate it. Thermodynamic projector completely eliminates this difficulty. From the 
other side, on the quasiequilibrium manifold with the moment parameterization (if it is 
found) it is easy to find the dynamics: simply write = J fiaJdq. The building of the 
thermodynamic projector may require some efforts. 

Finally, for each of the distributions \& it is easy to find its projection on the classical 
quasiequilibrium manifold \E' \E'^(.^^: it requires just calculation of the moments M(\E'). 
The analogue projection for the general thermodynamic projector is rather difficult: — > 
/ with the condition Pj'^(\E' — /) =0. This equation defines the projection of some 
neighborhood of the manifold Q on Q, but the solution of this equation is rather difficult. 
Fortunately, we need to build such operators only to analyze the fast processes of the 
initial relaxation layer, and it is not necessary to investigate the slow dynamics. 

9 A few words about the specifics of the computa- 
tional difficulties 

From the computational point of view, the main difficulties in realization of the described 
methods are in the calculation of the integrals of the form 



where are given functions of the vector q, ai, hi are the numbers, F is a function of 
one variable. The usual F are F{z) = e^]F{z) = 1/z,.... The usual dimension of Q in 
polymer physics is a few hundreds, number of different is a few dozens. 

In any case, the transition from the integration of the whole FPE to solution of the 
moment equations gives a considerable decrease of the computation time. 

In the methods of Legendre integrators and thermodynamic projector the computa- 
tional problems of the linear algebra are present: the solution of the system of linear 
equations Cfi = M (jHH), the problem of the orthogonalisation of vectors in Tf (jH!^ and 
so on. All these problems have the data which depends smooth on the current state of \E', 
and, consequently, on the time t. So, it is possible to solve these problems with the help 
of the perturbation theory and the methods of parametric continuation. These methods 
of the computational linear algebra are widely used and their details are well-known, so 
we are not discussing it here ( |2Sl I21j • 
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10 Accuracy estimation and postprocessing in invari- 
ant manifolds constructing 

Suppose that for the dynamical system the approximate invariant manifold has been 
constructed and the slow motion equations have been derived: 

dx I 

-^ = P,^,iJixsi)),XsieM, (65) 

where Px^i is the corresponding projector onto the tangent space Tj,.^, of M. Suppose that 
we have solved the system (fHHjl and have obtained Xsi{t). Let's consider the following two 
questions: 

• How well this solution approximates the real solution x{t) given the same initial 
conditions? 

• How is it possible to use the solution Xsi(t) for it's refinement without solving the 
system Q again? 

These two questions are interconnected. The first question states the problem of the 
accuracy estimation. The second one states the problem of postprocessing. 
The simplest ( "naive" ) estimation is given by the "invariance defect" : 

A., = (l-P,.JJ(a;,0 (66) 

compared with J{xsi). For example, this estimation is given by e = || Aa;^J|/|| J(xs/)|| using 
some appropriate norm. 

Probably, the most comprehensive answer to this question can be given by solving the 
following equation: 

^ = ^.^t) + DxJ{x)U^^^t)6x. (67) 

This linear equation describes the dynamics of the deviation 6x{t) = x{t) — Xsi{t) 
using the linear approximation. The solution with zero initial conditions 6x{0) = allows 
estimating Xsi robustness as well as the error value. Substituting Xsiit) for Xsi{t) + 6x{t) 
gives the required solution refinement. This dynamical postprocessing [2^] allows to refine 
the solution substantially and to estimate it's accuracy and robustness. However, the 
price for this is solving the equation ()67j) with variable coefficients. Thus, this dynamical 
postprocessing can be followed by a whole hierarchy of simplifications, both dynamical 
and static. Let's mention some of them, starting from the dynamical ones. 

1) Freezing the coefRcients. In the equation ()67|) the linear operator DxJ{x)\x^i{t) 
is replaced by it's value in some distinguished point x* (for example, in the equilibrium) 
or it is frozen somehow else. As a result, one gets the equation with constant coefficients 
and the explicit integration formula: 
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5x{t)= f exp{D*{t-r))A,^^^r^dT, (6^ 

where D* is the "frozen" operator and 6x{0) = 0. 

Another important way of freezing is substituting (|67|) for some model equation, i.e. 
substituting DxJ{x) for where r* is the relaxation time. In this case the formula for 
6x{t) has a very simple form: 

Sx{t) = [ e^A,^,(,)dr. (69) 



2) One-dimensional Galerkin-type approximation. Another "scalar" approxi- 
mation is given by projecting (jU7j) on A{t) = A^^j(t): 

AW. ^^l..<^^?^^. (TO) 

where ( |) is an appropriate scalar product which can depend on the point Xsi (for example, 
the entropic scalar product), D = DxJ{x)\x^i(t) or the self-adjoint linearization of this 
operator, or some approximation of it. 

The "hybrid" between equations (j7(J|) and (j67|) has the simplest form (but is more 
difficult for computation than eq. (fTOI) ): 

Here one uses the normalized matrix element instead of the linear operator D = 

D^J{x)\^^^(t)- 

Both equations (ffOjl and (ffTjl can be solved explicitly: 



5{t) = drexp (^j^ k{e)dej , (72) 
Sx{t) = A{T)dT exp (^j^ ki{e)de^ , (73) 

where m = hH) = i^. 

The projection of A^^^^ (t) on the slow motion is zero, hence, for post-processing analysis 
of the slow motion, the one-dimensional model ()7U|) should be supplemented by one more 
iteration: 

^^^^^ = 5(t)P.,,(,)(D.J(x,Kt)))(A(t)); 
Sx,i{t) = [ 6{T)P,^^^,^{D,J{xsiir))){A{r))dT. (74) 
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where 6(t) is the solution of (fTOj) . 

3) For a static postprocessing one uses stationary points of dynamical equations 
(IU7j) or their simphfied versions . (fTOj) . Instead of ^7\i one gets: 

DxJix)\^^^^t)Sx = -A^^^(t) (75) 

with one additional condition, Px^i^x = 0. This is exactly the iteration equation of the 
Newton's method in solving the invariance equation. 

The corresponding stationary problems for the model equations and for the projections 
of (j67|) on A are evident. We only mention that in the projection on A one gets a step of 
the relaxation method for the invariant manifold construction. 

For the static postprocessing with frozen parameters the "naive" estimation given by 
the "invariance defect" (j66|) makes sense. 

11 Example: Dumbell model, explosion of the Gaus- 
sian anzatz and polymer stretching in flow 

Here is an example of application of the thermodynamic projector method. In this ex- 
ample we consider the following simplest one-dimensional kinetic equation for the config- 
uration distribution function t), where q is the reduced vector connecting the beads 
of the dumbell. This equation is slightly different from the FPE considered above. It is 
nonlinear, because of the dependence of U on the moment M2[\l/] = / q^'^{q)dq. This de- 
pendence allows us to get the exact quasiequilibrium equations on M2, but this equations 
are not solving the problem: this quasiequilibrium manifold may become unstable when 
the flow is present j^H]- Here is this model: 

dt^ = -d,{a{t)q^>] + (76) 

Here 

a{t) = Kit) - ^/(M2(t)), (77) 

K{t) is the given time-independent velocity gradient, t is the reduced time, and the func- 
tion —fq is the reduced spring force. Function / may depend on the second moment of 
the distribution function M2 = J q^^{q,t)dq. In particular, the case / = 1 corresponds 
to the linear Hookean spring, while / = [1 — M2(t)/6]^^ corresponds to the self-consistent 
finite extension nonlinear elastic spring (the FENE-P model, first introduced in 150]). 
The second moment M2 occurs in the FENI5-P force / as the result of the pre-averaging 
approximation to the original FENE model (with nonlinear spring force / = [l — q^/h]~^). 
Leading to closed constitutive equations, the FENE-P model is frequently used in simu- 
lations of complex rheological flows as the reference for more sophisticated closures to the 
FENE model jHH 133 EHl- The parameter h changes the characteristics of the force law 
from Hookean at small extensions to a confining force for q^ — >■ h. Parameter h is roughly 
equal to the number of monomer units represented by the dumbell and should therefore 
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be a large number. In the limit 6 — * oo, the Hookean spring is recovered. Recently, 
it has been demonstrated that FENE-P model appears as first approximation within a 
systematic self-confident expansion of nonlinear forces [SH HE] • 

Equation (fTUI) describes an ensemble of non-interacting dumbells subject to a pseudo- 
elongational flow with flxed kinematics. As is well known, the Gaussian distribution 
function, 



exp 



2Mo 



solves equation (fTBj) provided the second moment M2 satisfles 



dM2 
dt 



1 + 2a{t)M2. 



(7J 



(79) 



Solution (fTSj) and (f?^ is the valid macroscopic description if all other solutions of the 
equation (fTUj) are rapidly attracted to the family of Gaussian distributions (fTHj) . In other 
words [3H], the special solution (fTSj) and (f?^ is the macroscopic description if equation 
(I78|) is the stable invariant manifold of the kinetic equation (j76j) . If not, then the Gaussian 
solution is just a member of the family of solutions, and equation (j79|) has no meaning of 
the macroscopic equation. Thus, the complete answer to the question of validity of the 
equation (f7^ as the macroscopic equation requires a study of dynamics in the neighbor- 
hood of the manifold (fTSj) . Because of the simplicity of the model (fTUj). this is possible to 
a satisfactory level even for M2-dependent spring forces. 

In the paper fI^ it was shown, that there is a possibility of "explosion" of the Gaussian 
manifold: with the small initial deviation from it, the solutions of the equation (fTUI) are 
very fast going far from, and then slowly come back to the stationary point which is 
located on the Gaussian manifold. The distribution function is stretched fast, but 
looses the Gaussian form, and after that the Gaussian form recovers slowly with the new 
value of M2. Let us describe briefly the results of 155]. 

Let = J q^""^dq denote the even moments (odd moments vanish by symmetry). 
We consider deviations fX2n = M2n - Mg, where = J q^^'^^dq are moments of the 
Gaussian distribution function (fTHj) . Let \E'(g, to) be the initial condition to the Eq. (f7S|l 
at time t = to- Introducing functions. 



P2n{t,to) = exp 



An I a{t')dt' 

to 



^0) 



where t > tg, and 2n > 2, the exact time evolution of the deviations /i2n ioi 2n> 2 reads 



and 



fJ'2n{t) 



lii{t) = P4(t,to)Ai4(to), 



fl2n{to) + 2n{An - 1) / fi2n~2{t')P2n{t' ,to)dt' 



P2n{t,tfi 



(81) 



(82) 



for 2n > 3. Equations ()80|). ()81|) and ()82j) describe evolution near the Gaussian solution 
for arbitrary initial condition \E'(g, to). Notice that explicit evaluation of the integral in 
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the Eq. ()80p requires solution to the moment equation ()79|) which is not available in the 
analytical form for the FENE-P model. 

It is straightforward to conclude that any solution with a non-Gaussian initial condi- 
tion converges to the Gaussian solution asymptotically as t ^ oo if 

lim / a{t')dt' < 0. (83) 

J to 

However, even if this asymptotic condition is met, deviations from the Gaussian solution 
may survive for considerable finite times. For example, if for some finite time T, the 
integral in the Eq. ()80|) is estimated as J^* a{t')dt' > a{t — t^), a > 0, t < T, then the 
Gaussian solution becomes exponentially unstable during this time interval. If this is the 
case, the moment equation (f7^ cannot be regarded as the macroscopic equation. Let us 
consider specific examples. 

For the Hookean spring (/ = 1) under a constant elongation (k = const), the Gaussian 
solution is exponentially stable for n < 0.5, and it becomes exponentially unstable for k, > 
0.5. The exponential instability in this case is accompanied by the well known breakdown 
of the solution to the Eq. (j79p due to infinite stretching of the dumbbell. Similar instability 
has been found numerically in three-dimensional flows for high Weissenberg numbers 

Eqs. (fTUjl and (|HT|) were integrated by the 5-th order Runge-Kutta method with adap- 
tive time step. The FENE-P parameter b was set equal to 50. The initial condition was 
0) = C(l — q^/hY^"^, where C is the normalization (the equilibrium of the FENE 
model, notoriously close to the FENE-P equilibrium 37J). For this initial condition, in 
particular, /i4(0) = —QlP'/[{h + 3)^(6 + 5)] which is about 4% of the value of M4 in the 
Gaussian equilibrium for h = 50. In Fig. ^ we demonstrate deviation /i4(t) as a function 
of time for several values of the flow. Function M2{t) is also given for comparison. For 
small enough k we find an adiabatic regime, that is fi^ relaxes exponentially to zero. For 
stronger flows, we observe an initial fast runaway from the invariant manifold with Ifi^l 
growing over three orders of magnitude compared to its initial value. After the maximum 
deviation has been reached, 114 relaxes to zero. This relaxation is exponential as soon as 
the solution to Eq. (f7^ approaches the steady state. However, the time constant for this 
exponential relaxation |aoo| is very small. Specifically, for large k, 

aoo = lima(t) = -^ + 0(ft-i). (84) 

Thus, the steady state solution is unique and Gaussian but the stronger is the flow, the 
larger is the initial runaway from the Gaussian solution, while the return to it thereafter 
becomes flow-independent. Our observation demonstrates that, though the stability con- 
dition ()83|1 is met, significant deviations from the Gaussian solution persist over the times 
when the solution of Eq. (f?^ is already reasonably close to the stationary state. If we 
accept the usually quoted physically reasonable minimal value of parameter b of the order 
20 then the minimal relaxation time is of order 40 in the reduced time units of Fig. ^ We 
should also stress that the two limits, k — > cx) and 6 — > 00, are not commutative, thus it is 
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Figure 1: Deviations of reduced moments from the Gaussian solution as a function of 
reduced time t in pseudo-elongation flow for the FENE-P model. Upper part: Reduced 
second moment X = M2/6. Lower part: Reduced deviation of fourth moment from 

1 /2 

Gaussian solution Y = — /i/ /b. Sohd: k = 2, dash-dot: k = 1, dash: k = 0.75, long 
dash: k = 0.5. (The figure from the paper [22], computed by P. Ilg.) 



not surprising that the estimation (j84j) does not reduce to the above mentioned Hookean 
result as 6 ^ 00. Finally, peculiarities of convergence to the Gaussian solution are even 
furthered if we consider more complicated (in particular, oscillating) flows K{t). 
In accordance with the anzatz for \E' can be suggested in the following form: 

= -4= (e-^+e-^). (85) 
Zay Zn \ / 

Natural inner coordinates on this manifold are a and Note, that now 7^ M2. The 
value 0"^ is a dispersion of one of the Gaussian summands in (jH3j) . 

To build the thermodynamic projector on the manifold ()85|) . the thermodynamic Lya- 
punov function is necessary. It is necessary to emphasize, that equations (fTB]) are non- 
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linear. For such equations, the arbitrarity in the choice of the thermodynamic Lyapunov 
function is much smaller. Nevertheless, such function exists. It is the free energy 



F = f/(M2[^]) -T5[^], 



^6) 



where 



^(In^ - l)dq, 



U(M2[^]) is the potential energy in the mean field approximation, T is the temperature 
(further we assume that T = 1). The thermodynamic properties of the mean field models 
in polymer physics are studied in the recent paper 

Note, that Kullback-form entropy S'fc = — J In {^) also has the form Sk = —F/T: 

^* = exp(-f/), 

Sk[^^J] = -{U)- j ^Ulii^ dq. 

If U{M2[^]) in the mean field approximation is the convex function of M2, then the free 
energy ()8fj|l is the convex functional too. 

For the FENE-P model U = - ln[l - M2/6]. 

In accordance to the thermodynamics the vector of flow of \1/ must be proportional to 
the gradient of the corresponding chemical potential /i: 



J= -5(^)V,/i, 
where /i = 4^, i? > 0. From the equation (jHHj) it follows, that 



^7) 



dU{M2j 2 



dMo 



+ In \1/ 



J 



If we suppose here i? = then we get 



J 



-D 



dU 
dM-y 



dm 
'dt 



1, 
2 
D 



q^ + -V,^ 



(89) 



When D = 1 this equations coincide with ()76|l in the absence of the flow: due to equation 
dHHl) dF/dt < 0. 

Let us construct the thermodynamic projector with the help of the thermodynamic 
Lyapunov function F (jHUj) . Corresponding entropic scalar product in the point has the 
form 

q f{q)dq ■ J q g{q)dq + J ^^^^ dq 



{f\9) 



dMi 



M2=M2[*] 



(90) 
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During the investigation of the anzatz (|H3j) the scalar product (|^. constructed for the 
corresponding point of the Gaussian manifold with M2 = cj^, will be used. It will let us 
to investigate the neighborhood of the Gaussian manifold (and to get all the results in 
the analytical form): 



M2=o-2 



q'fiq)dq ■ / q'giq)dq + av^ / f{q)g{q)dq (91) 



Also we will need to know the functional DF in the point of Gaussian manifold: 

1 



M2=(t2 



2(t2 



(92) 



(with the condition f f{q)dq = 0). The point 

dU{M2) 



dMo 



M2=<t' 



2 2a2' 



corresponds to the equilibrium. 

The tangent space to the manifold (|H3j) is spanned by the vectors 



4cTV27r 
1 

4^2^72^ 



5(^2) ' 

■ _(2±^(g + 0^-^x2 (g - ^)2 - a 

e 2<t2 ^ 2<t2 



21 



(93) 



_(2±o!g + ? (g-^)^ (g - c) 
-e 2<t2 1- e 2<t2 



0" 



cr 



The Gaussian entropy (free energy) production in the directions and 
simple form: 



DFAf,) = DFAU) 



dU{M2 



dM-y 



2(t2 



has a very 
(94) 



The linear subspace ker DF„2 in lin{ff^, f^} is spanned by the vector — /o-. 

Let us have the given vector field d'^ /dt = $(\1') in the point \E'({cr, ^}). We need to 
build the projection of $ onto the tangent space T„^^ in the point ^^({o", <;}): 



This equation means, that the equations for cr^ and vvill have the form 

da"^ dq"^ 
dt dt 



(95) 



(96) 
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Projection {(f^, ip^) can be found from the following two equations: 

^a + ^, = j q^^{^){q)dq- 

{^j^ + - h).^ = mnf^ - h).-. (97) 

where {f\g)„2 = (<l>(\l/)|/^ — /^)o-2, (pH) . First equation of (P7j) means, that the time 
derivative dM2/dt is the same for the initial and the reduced equations. Due to the 
formula for the dissipation of the free energy (jHH), this equality is equivalent to the 
persistence of the dissipation in the neighborhood of the Gaussian manifold. 

The second equation in ()97|1 means, that is projected orthogonally on ker DS To-,<^. 
Let us use the orthogonality with respect to the entropic scalar product (jHH). The solution 
of equations has the form 

da' ($1/. - A).. + M2($)((/,| - 



dt (/.-AI/.-A).2 

(98) 

_ -('^1/. - A).^ + M2(i')((M/.).2 - 

where $ = $(\l/), Mi{^) = J q'^{^)dq. 

It is easy to check, that the formulas are indeed defining the projector: if /„■ 
(or /^) is substituted there instead of the function $, then we will get ip^ = l,ip^ = Q 
(or (p„ = 0, 99^ = 1, respectively). Let us substitute the right part of the initial kinetic 
equations (fTBj) . calculated in the point \l/(g) = \1/({ct, <^}, g) (see the equation (|H3j)) in the 
equation instead of $. We will get the closed system of equations on cr^, in the 
neighborhood of the Gaussian manifold. 

This system describes the dynamics of the distribution function The distribution 
function is represented as the half-sum of two Gaussian distributions with the averages 
of distribution ±<^ and mean-square deviations a. All integrals in the right-hand part of 
are possible to calculate analytically. 

Basis (/cr, /^) is convenient to use everywhere, except the points in the Gaussian man- 
ifold, ^ = 0, because if <^ — > 0, then 

To analyze the relaxation in the small neighborhood of the Gaussian manifold it is more 
convenient to use another basis: 



It corresponds to a reparametrization of the initial manifold (jH3j) : 



,2 

;^(/a' "~ A)- 



^ ^^^/|^ [e-^(^^ + e-5(e^ ) . (99) 
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Let us analyze the stability of the Gaussian manifold to the "dissociation" of the Gaussian 
peak in two peaks ()85|). To do this, it is necessary to find first nonzero term in the Taylor 
expansion in of the right-hand side of the second equation in the system ()98p. The 
denominator has the order of <^^, the numerator has, as it is easy to see, the order not less, 
than (because the Gaussian manifold is invariant with respect to the initial system). 



Let us denote Ga 



e 



Then we get 
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Let us calculate = $(^({cr, with the accuracy up to 

1 



M2(^a|vl/({a,,})) = l; 



-a9g(g*({or,^})) = aG<,(g) 



2a4 



(7^ 



a 



4 U4 a6 



n g^ 7g^ 
24^^ 8^ 



+ 



M2(-a9,(gvI/({a, ^}))) = 2a(or2 + + o I 



The diffusion part gives the zero contribution to the numerator of the equation (|98|): 

-(AIA-A) + (AIA-A)=o, 

therefore to find dq/dt it is sufficient to use $1 = — «9q(g\E'), so we get 
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4 

4 
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Thus 



^-rr = + o — r 
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(100) 
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Figure 2: Phase trajectories for two-peak approximation, FENE-P model. The vertical 
axis (<j = 0) corresponds to the Gaussian manifold. The triangle with a{M2) > is the 
domain of exponential instability. 

So, if a > 0, then grows exponentially ~ e"*) and the Gaussian manifold is unstable; 
if a < 0, then decreases exponentially and the Gaussian manifold is stable. 

The form of the phase trajectories is shown qualitative on the figure El 

For the real equation FPE (for example, with the FENE potential) the motion in 
presence of the flow can be represented as the motion in the effective potential well 
U{q) = U{q) — Kq"^. Different variants of the phase portrait for the FENE potential are 
present on the figure El Instability and dissociation of the unimodal distribution functions 
("peaks") for the FPE is the general effect when the flow is present. 

The instability occurs when the matrix d'^U/dqidqj starts to have negative eigenvalues 
{U is the effective potential energy, U{q) = U{q) — ^ij'ii'ij)- 

The stationary polymodal distribution corresponds to the persistence of several local 
minima of the function U{q). The multidimensional case is different from one-dimensional 
because it has the huge amount of possible configurations. All normal forms of the 
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Figure 3: Phase trajectories for two-peak approximation, FENE model: a) A stable 
equilibrium on the vertical axis, one stable peak; b) A stable equilibrium with ^ > 0, 
stable two-peak configuration. 



catastrophe of "birth of the critical point" are well investigated and known ^3]. Every 
dissociation of the peak is connected with such a catastrophe. The number of the new 
peaks is equal to the number of the new local minima of U. 

It is not very difficult to perform the analysis of the equations for every quantity 
of peaks and every potential. Moreover, for the polynomial potentials all the necessary 
integrals are possible to calculate analytically (if the coefficients of the scalar product and 
entropy production are taken in the Gaussian point). The same situation is also for the 
general Gaussian distributions: 

Cs = const ■ exp ^-^ ^(S"^)^^?^^^^ , 

where S is the covariance matrix. Here in the equation for the effective energy we have 
the symmetric part of the tensor = d'^U / dqidgj . The presence of the unsimmetric part 
may lead to the relaxation oscillations (both for FPE and for the peak dynamics). 

For the modeling of dynamics of the multimodal distributions for FPE with the pres- 
ence of the flow (the flow may be nonstationary) it seems to be useful to use the physically 
clear modeling of the distribution function as a sum of the flnite number of the Gaussian 
peaks. Thermodynamic projector gives us an opportunity to make this models thermo- 
dynamically consistent. 
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12 Conclusion 



In this work we presented a toolbox for the development and reduction of the dynamical 
models of nonequilibrium systems with the persistence of the correct dissipation. 

The basic notions of this toolbox are: entropy, quasiequilibrium (MaxEnt) distribution, 
dual variables, thermodynamic projector. 

The main technical ideas are: Legendre Integrators, dynamical postprocessing, trans- 
formation of almost arbitrary anzatz to a thermodynamically consistent model via ther- 
modynamic projector. 

The Legendre Integrators are based on a simple, but very useful idea: to write and 
solve dynamic equations for dual variables. This idea is efficient, because to obtain the 
dynamic equations for dual variables it is necessary to solve linear equations. To get 
the usual quasiequilibrium dynamical equations for the moments, we should solve nonlin- 
ear (transcendent) equations. Sometimes it happens that these equations can be written 
down in the explicit form (Vlasov equation, Euler equation, ten moments Gaussian ap- 
proximation in gas kinetics [El ESI); t>ut usually these equations remain in implicit form 
with right-hand sides derived by a system of transcendent equations. 

The post-processing is necessary for accuracy estimation. It gives us the cheapest way 
to improve the solution obtained by the Legendre Integrators. 

Termo dynamic projector allows to transform almost arbitrary anzatz into a physically 
consistent dynamic model with persistence of dissipation. The simplest example, discussed 
in details, is the two peaks model for Gaussian manifold instability in polymer dynamics. 
This type of models opens a way to create the computational models for the "molecular 
individualism" HQ EH HE] • 

The simplest model of the molecular individualism is the "Gaussian parallelepiped". 
The distribution function is represented as a sum of 2"^ Gaussian peaks located in the 
vertixes of centrally symmetrical parallelepiped: 



^ (4 • • 

where n is dimension of configuration space, 2<^j is the vector of the ith edge of the 
parallelepiped, S is the one peak covariance matrix (in this model E is the same for all 
peaks). The macroscopic variables for this model are: 

1. The covariance matrix S; 

2. The set of vectors q (or the parallelepiped edges). 

The dimension is n{n + l)/2 + mn. 

The number m {m < n) is the estimated number of nonstable directions of motion 
(dimension of instability). To include the nongaussian equilibrium the "Gaussian paral- 
lelepiped" should be deformed to nongaussian "peaks parallelepiped". Technical details 
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will be discussed in the separate paper. The structure of "peaks parallelepiped" leads to 
the molecular individualism in such a way: each individual molecule belongs to a domain 
of a peak in configuration space. The number of these peaks grows significantly with the 
dimension of instability, as 2"*, and even if m = 3, than the number of peaks is 8, and 
one should discover 8 distinguished sorts of molecular configurations. On the other hand, 
in projection on a line this amount of peaks can form a distribution without a sign about 
peak structure, hence, the study of properties of ensembles (viscosity, stress coefficient, 
etc.) can be without any hint to a cluster structure in configuration space. 
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